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We present a very brief description of the Hartree-Fock method in nuclear structure 
physics, discuss the numerical methods used to solve the self-consistent equations, 
and analyze the precision and convergence properties of solutions. As an appli- 
cation we present results pertaining to quadrupole moments and single-particle 
quadrupole polarizations in superdeformed nuclei with A~60. 



1 Introduction 

Self-consistent methods have been used in the low-energy nuclear structure 
studies over many AfparSj-and represent a mature field with numerous suc- 
cessful applicationsEmaQO A number of computer codes solving the nuclear 
Hartree-Fock (HF) problem have already been developed. Two types of effec- 
tive nucleon-nucleon interactions have been mainly employed. Starting with 
the work of Vautherin and Brink □ many authors have applied the nucle 
HF theory with the Skyrme effective interaction, while the work of Gogny t 
initiated numerous studies with the force which carries his name. 

The methods employed to solve the HF equations depend mainly on the 
effective force used and on the assumed symmetries of the many-body wave 
functions. For the solutions which allow at least triaxial deformations, two 
different methods have been applied for the two above mentioned effective in- 
teractions. The first one, used in conjuncture with the Skyrme interaction, is 
formulatecLin the spatial coordi nates and makes use of the finite-differencep 
or Fourierl^J or spline-collocationEJI13 methods to approximate differential op- 
erators. Ihe solution is then obtained by using the imaginary time evolution 
operatorE3 

The second one, used for the finite-.range Gogny interaction, employs a 
truncated harmonic oscillator (HO) basisotJ and solves the problem either by. 
an iterative diagonalization of the mean-field Hamiltonian, or by the gradientliS 
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or the conjugate gradient0 methods. Recently, the method which incorporates 
the advantages of both existing approaches, and combines the robustness of 
the Cartesian_HO basis with the simplicity of the Skyrme interaction, has been 
implemented E£l 

The methods using spatial coordinates have several advantages. First of 
all, various nuclear shapes can be easily treated on the same footing; the same 
cubic lattice of points in three spatial dimensions is suitable to accommodate 
wave functions with, in principle, arbitrary deformations restricted only by 
a specific symmetrization of the lattice. This allows easy studies of systems 
for which the deformation is not a priori known, or is ill defined because 
of deformation instabilities or a shape coexistence. Secondly, using spatial 
coordinates allows studies which address the asymptotic form of nucleonic 
wave functions at large distances. This is particularly important for a precise 
description of weakly bound nuclei, where the use of spatial coordinates is a 
necessitytS Third, for the Skyrme zero-range interaction, the mean fields are 
local (apart from a velocity dependence) and can be easily programmed in the 
spatial coordinates. Last but not least, the treatment of wave functions on large 
lattices (12x12x12 is a typical example) is easily amenable to vectorization or 
parallelization of the algorithm. 

Methods using the HO basis have other advantages. Firstly, the basis 
provides a natural cut-off for many operators which otherwise are unbound 
and require particularly delicate treatment in the spatial coordinates. This 
concerns in particular the multipole moment and the angular momentum op- 
erators which are often used as constraining operators. For the corresponding 
constraints the solutions can become unstable when the non-zero probability 
amplitudes (wave functions) move towards large distances as it is the case 
for e.g. weakly bound nucleons. Secondly, much smaller spaces are usually 
required to describe the nuclear wave functions within a given precision. Typi- 
cally a basis of about 300 HO wave functions is sufficient for most applications. 
Third, the iterative diagonalization of the mean field Hamiltonian can be used 
to find the self-consistent solutions, which provides a rapidly converging al- 
gorithm, and, last but not least, scalar or superscalar computers can also be 
used, because the typical sizes of the information handled is smaller and the 
performance is less dependent on the use of a vector processor. 

Detailed discussion of the stability, convergence, and efficiency of methods 



s recently been presented 
In the present communi- 



using the expansion on the Cartesian HO basis ' 
together with the description of the HFODD code 
cation we give a brief resume of the methods employed (Sec.0), present tests 
of the precision and discuss the CPU times required (Sec. Q), and illustrate 
possible applications by presenting a few recent results obtained for the A~60 
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nuclei (Sec. ||). 



2 Hartree-Fock method 

The eigenequations for the HF single-particle Routhians h' a are called the HF 
equations, 

h' a ip it(X (ra) = e' iia i/>i, a (ra), (1) 
where the neutron (a=n) and proton (a=p) Routhians read 

- r 2 



//,', = -|-A + T a + 6 ap U COUl + C/ mUlt - IjJyJy (2) 



Detailed expressions for the nuclear mean field operators T a , Coulomh_poten- 
tial U Conl , and multipole constraint potential [/ mult are given in Refta After 
solving Eq. ([!]) one calculates the total energy of the system as a sum of the 
kinetic, Skyrme, Coulomb, and pairing terms 

£ = £ kin + £ SkyTme + £ Coul + £ pair . (3) 

The same energy can be calculated by using the Routhian eigenvalues e[ and 
occupation probabilities vf a as 

c 1 \ ^ „.2 I _|_ 1_ ckin , cpair crear , 1 cCoul cmult ccran f a\ 

c ~ 2 / j u i,a L i,a ' 2 ' 3 C cxch c corr c corr ■ \^) 



Here, the rearrangement energy £ rcar results from the density dependence of 
the Skyrme interaction. Similarly, the Coulomb exchange energy £^°^ can be 
considered as resulting from a zero-order interaction term depending on the 
density as p p ' . The remaining terms correct for the fact that Routhians 
contain constraint terms which are not present in the total energy. 

The difference between the energies £ and £ is exactly equal to zero when 
the densities and fields do not change from one iteration to the next one. Hence 
their difference 

5£ = £-£ (5) 

provides a useful measure of the quality of convergence, and is called the sta- 
bility of the HF energy. 

The HF equations (|) can be solved by expanding the single-particle wave 
functions ipi(ra) onto the deformed HO wave functions ^> 7 i, J .n y n^.s,{ riJ ) m the 
Cartesian coordinates, i.e., 

N x N v AT, 

^h=EEE E ^ w *vw*,s*(™)' (6) 

n x =0 n y =0n z =0 Sz= - 1 ,| 
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Here N Xl N y , and N z are the maximum numbers of the HO quanta corre- 
sponding to the three Cartesian directions. However, the sums over n x , n y , 
and n z are performed over the grid of points which form a pyramid rather than 
a cube. 

The HO wave functions have the standard form 

^Pn x nyn z ,s z 

(ra) = ip n:c (x)ip nv (y)^n z {z)S Sz(T , (7) 

where 

^{x») = blH£}{^)e-^l, (8) 
and ^fi—b^Xfj, are dimensionless variables scaled by the oscillator constants 



ftp = ^mu^/h. (9) 

Polynomials (£) are proportional to the standard Hermite orthogonal poly- 
nomials ff„(£)@ 

Hji 0) (O = (V5F2"n!)-^ n (O. (10) 

Accuracy of the solution of the HF equations with the wave functions ex- 
panded onto the Cartesian HO basis, Eq. (0), depends on the three parameters 
hu) x , Tiujy, hoj z defining the HO frequencies in three Cartesian directions, and 
on the number M of the HO states. included in the basis. In the code hfodd 
we use the standard prescription EiiE3 to chose the HO states included in the 
basis, namely, the M states with the lowest HO single-particle energies, 

e^rijn, = huj x (n x + |) + huj y (n y + |) + hu> z (n z + |), (11) 

are selected among those which have u x <Nq, n y <No, and n z <No, where iVo is 
the fixed maximum number of HO quanta. It should be noted that in general 
both M and No have to be specified to define the basis. Only for large No, 
the basis is defined solely by M and does not depend on iVo- In this case, 
the grid of points (n x , n y ,n z ) defining the states included in the basis forms a 
pyramid in three dimensions, with the inclined face delimited by the condition 
£n x n n,<const. On the other hand, only for small values of No the basis is 
defined solely by No and does not depend on the energy cut-off. In this case 
the corresponding grid of points n x n y n z forms a cube of the size No- In all 
intermediate cases the shape of the basis corresponds to a pyramid with the 
corners cut off, or to a cube with the corners cut off. Usually No is chosen 
large enough so that all the states allowed by the energy cut-off are included 
in the basis. 
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3 Precision and convergence 



The HF equations (|lj) can be solved by using standard iterative methodsld which 
are based on the self-consistency principle. Namely, the self-consistent solution 
is attained when the mean fields T a appearing in the mean-field Routhian h' a 
equal to the mean fields calculated from the occupied single-particle states 
Tpi,a{r<j). Therefore, by iteratively calculating the mean fields and using them 
to obtain new single-particle states one may eventually reach the self-consistent 
solution. A convergence of such a procedure is not guaranteed and, in fact, 
in nuclear applications the method such as formulated above almost always 
diverges. The reason for that is the fact that the consecutive corrections of the 
mean fields turn out to be too violent and lead to too strong modifications of 
the states. Since early years of practical implementations of the HF method 
in nuclear structure physics the simple method to remedy this situation has 
been devised. Namely, only the fraction / of the mean fields calculated from 
the single-particle states is used in every iteration, and it is combined with the 
fraction 1—/ of the old mean fields. 

In practical calculations, the optimal value of / is determined by checking 
the convergence rate. An example of such an analysis is presented in Fig. |], 
where the stability energy (^) is shown for the values of / ranging from 30% 
to 90%. It can be seen that the rate of convergence increases very fast when 
/ increases from 30% to about 80%, but then it decreases again dramatically. 
Unfortunately, the optimal value of / depends critically on the shell structure 
of the nucleus; it is rather high for magic nuclei (as superdeformed ggDyge) 
which have large shell gaps, and has to be lowered when there are several 
single-particle states near the Fermi surface. In practical calculations it is more 
efficient to use rather low values of /, and ensure convergence independently of 
the shell structure, than to risk a divergence. Usually, f—50% can be a safely 
recommended value. 

The total CPU times required to solve the HF equations depend not only 
on the convergence rate governed by /, but also (of course) on the overall 
speed of the computer and (very importantly) on the performances of the spe- 
cific compiler optimization procedures. This is illustrated in Fig. ||, where the 
relative times required to perform specific tasks are shown in percent of the 
total time. The total CPU times required to calculate the complete superdefor- 
med band in ^Dyse range from one hour on the Cray C-90 supercomputer to 
several hours on fast workstations (SG Power Challenge L or IBM RS/6000), 
and to about 36 hours on a typical workstation (SUN-Hyper). However, these 
times result from averaging very different performances pertaining to specific 
tasks. Typically, three tasks take most of the time, namely, the calculation of 
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Figure 1: Stability energy as function of number of iterations for the calculations performed 
for the superdeformed state of ^gDyse- Different curves correspond to different values of 
the fraction / (see text), and lead to different numbers of iterations which are necessary to 
improve the stability by an order of magnitude. 



the HO matrix elements of the mean fields, the diagonalization of the Routhi- 
ans, and the calculation of the Coulomb field. On a vector computer, the first 
of these tasks takes, relatively long time (40%) because it is composed of mul- 
tiple short loops, LL3 whereas the last one contributes very little (5%) due to 
rather long loops it requires. On superscalar computers these three tasks take 
almost the same times of about 20% each. The example of typical workstation 
shown in Fig. || illustrates the importance of the optimization options; here 
the compiler failed to optimize the subroutine calculating the mean fields and 
therefore this task takes now almost 30% of the CPU time. 



4 Superdeformed bands in the A^60 nuclei 

Very recently there have been several reports p f superdeformed rotational 
bands discovered in nuclei around A~60EjO't3r3 Specific nuclei in this region 
have already been addressed fcom the point of view of the Nilsson-Strutinskyc3 
and Relativistic Mean Field c3 descriptions. The single-particle structure of 
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Figure 2: Percentages of the CPU time spent on performing separate tasks required by the 
HFODD code to solve the HF equations. 



the light superdeformed nuclei is very simple and thus amenable to a system- 
atic description in terms of a microscopic self-consistent mean-field approach. 
Using the code HFODD (vl.25)p hundreds of bands have recently been calcu- 
lated in light Ni-Ga nuclei. c3 Based on these results one can find and analyze 
the generic features exhibited by rapidly rotating nuclei. In the present com- 
munication we report on quadrupole polarizabilities characterizing individual 
single-particle orbitals in this region. ._. 

Following the method introduced in the A^150 nuclei we aim at a de- 
scription of the proton quadrupole moments Q20 in terms of contributions 
given by the particle (p) and hole (h) states defined with respect to the magic 
superdeformed nucleus |oZn3o, i.e., 

Q20 = Q20 [loZnao] + 6 ®20 n p + E (12) 

p h 
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Figure 3: Deviations of calculated HF proton q uad rupole moments in A~60 nuclei from the 
sum of individual single-particle contributions fll2|). 



The core gpZi^o corresponds to all n p — and all rih — 0, while the other superde- 
formed states in this region can be obtained by creating up to 6 particles with 
rt p =l and/or up to 6 holes with nh=l, both for protons and for neutrons. Ej 
Hence, the sum in Eq. (|l2|) contains altogether 24 unknown coefficients SQ20 
which are fitted by the least-square method to reproduce all calculated proton 
quadrupole moments Q20J separately at every rotational frequency huj. 

The quality of describing the proton quadrupole moments by the simple 
additive contributions is illustrated in Fig. || where we show deviations of the 
calculated HF values from the sum given in Eq. (H). One can see that for most 
of the bands the additivity hypothesis is correct up to 0.05 b. This constitutes 
a remarkable agreement with the values of Q20 which are of the order of 2-3 b. 

In Fig. ^ the values of individual contributions SQ20 are shown as functions 
of the rotational frequency for the 24 orbitals considered. It seen that both 
proton and neutron orbitals significantly contribute to the proton quadrupole 
moment. This is so because the direct contributions of protons are of the sim- 
ilar order as the polarization contributions existing for both kinds of particles. 
The positive-parity orbitals [440] 1/2 and [431]3/2, originating from the lgg/2 
intruder orbital, have a larger influence on the quadrupole moment than the 
negative-parity orbitals (note that the scale is twice expanded in the top panel 
of Fig. Deformations gradually increase when consecutive 1 179/2 orbitals 
are occupied. The influence of proton positive-parity orbitals is about twice 
larger than that of their neutron counterparts. The negative-parity [303] 7/2 
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Figure 4: Contributions <5Q20 of individual single-particle orbitals to the proton quadrupolc 
moments of nuclei in the A~60 region. The orbitals are denoted by the dominant Nilsson 
labels. Solid (dashed) lines correspond to particle (hole) states, while the closed and open 
symbols correspond to positive (r=i) and negative (r=— i) signatures, respectively. Error 
bars are determined by the least-square procedure. 
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orbitals, originating from the I/7/2 extruder orbital, significantly contribute to 
the quadrupole moment too. In this case the neutron and proton contributions 
are almost equal, i.e., the direct proton contribution is very small. On the other 
hand, the proton contribution from the [312]5/2 orbital is much larger than 
that of the neutron counterpart; hence this orbital contributes mostly through 
the direct quadrupole moment. 

Contributions coming from the signature-partner orbitals are almost equal, 
except from those related to the [431]3/2 and [310]l/2 orbitals where the 
positive-signature partners contribute significantly more. These two orbitals, 
and also [303] 7/2, exhibit contributions slightly dependent on the angular fre- 
quency. However, these variations are relatively much smaller (do not exceed 
0.2 b) than the variation of the core quadrupole moment Q20 [30^30]. In- 
deed, the latter decreases quite substantially with increasing Jfcj, from 3.28 b 
at 7iw=0.2MeV to 2.56 b at 7iu;=1.6MeV. This is at variance with the results 
calculated for efDy86> where the core quadrupole moment is fairly indepen- 
dent of spin. 
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